function [newP origIndex newIndex] = interpolateP(P, sites)
% linear interpolation proportional to length  
disp('Linear Interpolation of Sites....');
newP(1,:) = P(1,:);
origIndex = [1];

for i=2:size(P,1)
    distance = norm(P(i-1,:) - P(i,:));
    M = floor(distance/2000);
  
    for j=1:M-1
       newP = [newP; ((M-j)*P(i-1,:)+j*P(i,:))/M]; 
    end
    newP = [newP; P(i,:)];
  
    origIndex = [origIndex; size(newP,1)];
end

tempIndex = 1:size(newP,1);
for i=1:length(origIndex)
    tempIndex(origIndex(i)) = 0;
end
newIndex = tempIndex(tempIndex > 0);

% adjust height of interpolated sites according to terrain
disp('Adjusting height of interpolated sites');

for i=1:length(newIndex)
    dist = [];
    for j=1:size(sites,1)
        distVertex(j) = norm(newP(newIndex(i),[1 2]) - sites(j,[1 2]),2); 
    end

    % among closest 10 vertices pick one with lowest height
    [distVal distIndex] = sort(distVertex);
    closestVertices = sites(distIndex(1:10),:);
    
    [heightVal heightIndex] = sort(closestVertices(:,3));
    
    newP(newIndex(i),:) = closestVertices(min(heightIndex),:);

end
